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IN THE UNITED STATES PATENT AND TRADEMARK OFFICE 



APPLICATION FOR LETTERS PATENT 



it it ic 



METHOD FOR PREDICTING DYNAMIC 
PARAMETERS OF FLUIDS 
IN A SUBTERRANEAN RESERVOIR 



INVENTORS 

VLADIMIR B. PISETSKI 
VALERI V. KORMILCEV 
and 

ALEKSANDER N. RATUSHNAK 



ATTORNEY'S DOCKET NO. TR01-P04 



1 METHOD FOR PREDICTING DYNAMIC PARAMETERS OF FLUIDS 

2 IN A SUBTERRANEAN RESERVOIR 

3 CROSS REFERENCE TO RELATED APPLICATION 

4 This application is a Continuation in Part of U.S. Patent Application Serial Number 

5 09/134,616, filed August 14, 1998, which is a Continuation in Part of U.S. Patent Application 

6 Serial Number 08/909,454, filed August 11, 1997, now U.S. Patent No. 5,796,678, issued August 

7 18, 1998. 
8 

9 FIELD OF INVENTION 

10 This invention pertains to a method for predicting the dynamic parameters of fluids in a 

a 1 subterranean reservoir. 

%3 BACKGROUND OF THE INVENTION 

." "l4 In the search for subterranean fluids (typically natural gas and/or petroleum), classical 

li5 methods have used seismic imaging to generate an image of a subterranean region of exploration. 

^16 Modem three dimensional (3D) seismic imaging has aided a great deal in developing credible 

H? models of the region of interest, hiformation regarding travel time of seismic signals can be used to 

18 estimate velocities of various subsurface layers, which provides information indicative of the rock 

19 type. Hydrocarbons are often trapped by geologic faults, and so the seismic images are 

20 subsequently interpreted to identify rock types which, in conjunction v^th geologic faults or other 

21 phenomenon, would result in an accumulation of fluids in an area in the region of study. However 

22 seismic information provides only limited information. Unless the geologic fault has a seal, fluids 

23 will not be trapped. Unfortunately, seismic survey data does not contain sufficient information to 

24 indicate the presence or absence of a seal. Other information is desirable which, if used in 
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1 conjunction with seismic images, would improve the probability of finding hydrocarbons. 

2 It is known that the field of so-called prestresses in the earth's crust is determined by two 

3 systems of external forces, namely, gravitational and tectonic forces. The joint action of these forces 

4 can result in compaction as well as in unloading of rock masses with respect to their normal 

5 gravitational or lithostatic stress condition. The existing theory describing stress fields in the earth's 

6 crust is absolutely correct for elastic and continuous media. At the same time information derived 

7 fi-om a number of wells, drilled down to 10 km, points to the presence of fractures in sedimentary 

8 and crystalline parts of the earth's crust. Since there is little reason to expect presence of spatially 

9 localized fi^actures, we can affirm that any real rock unit is bounded by a closed system of firactures 
Mo and represents a discrete system. Taking into account a well-known concept that any basin 
^ represents a tectonically young and consequently active structure, these are all reasons to consider 
ll2 any imit of the sedimentary cover as a part of discrete dynamic system. To date, no one has thought 
uB to use information pertaining to stress in the earth's crust as part of a method for determining the 
H4 presence of fluid in a subterranean formation. Our invention uses seismic imaging information in 

5 conjunction with formation stress information to improve the reliability of seismic exploration for 

%6 hydrocarbons. 

17 The pubUcation of J. D. Byerlee ("Friction of Rocks", Pure Applied Geophysics, 116 

18 (1978), pp. 615-626) is usefiil for estimating parameters of stress condition of jfractured and discrete 

19 rock units. He has found that the maximum differential stress (the difference of main stresses) in the 

20 upper part of the earth's crust is limited by the shear strength of firactured rocks. Laboratory 

21 investigations show that the sliding fiiction on the surface of fi:'actured rocks exists until the external 

22 load reaches the critical value of brittle failure. At the same time the strength of such a discrete 

23 rock unit is determined by cohesion of fractures on its surface and does not practically depend on 

24 the elastic moduli of the rock, temperature, strain value, and the type of sliding surface. It is 
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1 determined that the cohesive force r is a linear function of normalload (jn (Byerlee's law): 

2 r = 0.85 cTn 3 < o-n < 200 Mna (1) 

3 or in the terms of main stresses, the limit value of horizontal stresses down to the depths of 5 km 

4 is estimated to be 

5 cri«5cr3 o-3<110Mna (2) 

6 In other words the main component of stress condition variation in discrete media is 

7 associated with the change of vertical load. Lisignificant movement in the base of a sedimentary 

8 basin (for example, in the crystalline basement) results in considerable change of the horizontal 
^ 9 component of stress. The strength limit of a rock unit is rapidly reached, and all the sedimentary 
'HO cover is set in motion continuously approaching the isostatic equilibrium. Thus, it is reasonable to 
^il suggest that real sedimentary formations have a high degree of mobility in vertical direction. As 
1^2 time of stress relaxation characteristic for discrete media is significant, the stress distribution in the 
- 13 basin has distinctive contrasting character and is governed by the principles of block dynamics. 

ry4 The most fundamental theoretical results in the field of elastic wave propagation in 

□5 prestressed media were derived by M. A. Biot (textbook, Mechanics of Incremental Deformations, 

He New York, 1965). Biot introduced a new type of strain, strain of soUd rotation. Its introduction is 

17 related to non-hydrostatic character of prestress field in the medium (i.e. horizontal and vertical 

18 components of the stress field are not equal to each other). In this case each unit volume of the 

19 medium has been rotated through some angle with respect to its initial position in the unstressed 

20 medium. As a result additional strains appear in the wave field, and their amplitudes are directly 

21 proportional to the difference between horizontal and vertical components of the prestress field {ai 

22 ^ a 3)' A number of researchers studied the equations derived by M. A. Biot and came to the 

23 conclusion that the influence of initial stresses on traveltime and, especially, amplitude parameters 
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1 of seismic waves^ can be very significant, and the degree of this influence is proportional to the 

2 ratio P/^ (where P is pressure and ^ is the shear modulus). 

3 A significant number of publications are related to investigation of the stress condition of 

4 the earth with the aid of seismic methods. But practically all these investigations were based on the 

5 analysis of traveltime parameters and were reduced to attempts to explain complex distribution of 

6 velocity parameters by the influence of lithostatic pressure. In this case the stress condition of the 

7 subsurface is taken into account indirectly using theoretically and experimentally derived 

8 relationships between velocity and external load. In low-pressure regions the amplitudes and 

9 frequencies of the reflected signals completely depend on the density of fractures and the load 
irtO applied to them. And the earth material (its elastic moduli) significantly influences the amplitude 
f h 1 parameters of the signals only at the loads close to the strength limit of discrete media. 

Hi 2 What is needed then is a method for improving the ability to discover hydrocarbons in a 

'Ms subterranean formation and which makes use of available data, an in particular seismic data which 

\i\4 is likely to be collected as part of a exploratory effort in any event. 
Us 

m SUMMARY OF THE INVENTION 

17 The presented invention has been developed as a method for seismic data interpretation 

18 targeted at discovery of the areas of the accumulations of the fluids. We have discovered a 

19 relationship between reflection coefficient and general rock pressure, which is the basis of our 

20 invention. It is assumed that the real sedimentary sequence is a discrete medium subjected to 

2 1 inhomogeneous stresses created by two types of the forces, namely gravitational and tectonic forces. 

22 As a result of a sum influence of these two forces, the basin, as a discrete system, is found in 

23 continuous movement or in the condition of inhomogeneous stresses. In correspondence with this 
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1 point of view any stratigraphic sequence is found in different stress conditions in lateral direction, 

2 and consequently the amplitude and frequency parameters of reflected signals depend on the value 

3 of rock pressure acting on given element of a seismic reflector. 

4 The first step of the developed interpretation method consists in tracking one or several 

5 target reflections along stacked time sections. A standard tracking procedure and standard 2D or 3D 

6 processing flov^s are used. The pressure gradients are found using picked travel times for tracked 

7 seismic horizons and by calculating parameters of Hillbert transformation (instantaneous 

8 amplitudes and phases) for each identified signal. The interpretation formula is directed to 

9 estimation of relative pressure changes under the assumption that the smoothed values of 
instantaneous amplitudes and derivatives of instantaneous phases (instantaneous frequency) 

r|l correspond to general normal or gravitational pressure within all the study area, while all deviations 

fl2 from the average value represent relative estimates of anomalous pressures and correspond to local 

ffl3 regions of decompression and compression. 

H;4 The obtained map of relative pressure gradients is laid over the isochron map of a tracked 

5 horizon improved using the instantaneous phase picking. The principle of determining traveltimes 

^36 more accurately consists in correction for signal frequency variation along the reflecting horizon in 

17 correspondence with the derived dependence of the frequency from pressure. The resulting map of 

18 relative pressure variation in combination with reflection time map represents the basis for 

19 identifying most probable locations of fluid accumulations. Localized region of anomalous 

20 unloading coincident with a closed or semi-closed time high is a physically valid area of 

2 1 hydrocarbon and water accumulation. 

22 The information from the resultmg map of relative pressure variation in combination with 

23 reflection time map can be further combined with classical seismic interpretation techniques to 

24 further enhance the prediction of the presence of fluids in a subterranean formation. 
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1 While the invention described herein is particularly useful for predicting the presence of 

2 fluid accumulations in a subterranean reservoir, the relationship between reflection coefficient and 

3 general rock pressure can also be used to predict rock pressure or pressure gradients within a 

4 subterranean formation. Such information can be useful in and of itself, for example, in the study 

5 and prediction of earthquakes. 

6 The invention further includes a method for predicting the dynamic parameters of fluids in 

7 a subterranean reservoir. More specifically, the dynamic parameters of such fluids are 

8 determined by determining the relative estimate of the total ground pressure on the attributes of 

9 seismic signals related to a reflecting horizon, or any selected interval of the seismic section, 
yo This provides for a method for determining the location of the accumulation fluids in a 
Ml subterranean formation. The method includes the steps of determining a first velocity vector 
M2 "Vx" for migration of fluid in a region of interest in the subterranean formation. The first 
ffl3 velocity vector includes attributes of speed and direction of flow of fluid in a first direction in the 
^^4 region of interest. The method further includes determining a second velocity vector "Vy" for 
J j5 migration of fluid in the region of interest. The second velocity vector includes attributes of 
5 6 speed and direction of flow of fluid in a second direction in the region of interest. The velocity 

17 vectors are then extrapolated to identify the fluid accumulation location. The first and second 

18 velocity vectors are primarily functions of supplementary pressure "dP" in the region of interest, 

19 the permeability "c" of the region of interest, and the viscosity "u" of the fluid in the region of 

20 interest. 

21 The supplementary pressure can be determined by identifying pressure gradients within the 

22 region, the region being characterized by a seismic image of a stacked time section representing 

23 horizons within the region. The permeability of the media within the region, and the viscosity of 

24 the fluid within the region, can either be determined mathematically or fiom geological data. 
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1 In addition to the methods disclosed and described herein, our invention further includes 

2 computer apparatus for implementing the methods. 
3 

4 BRIEF DESCRIPTION OF THE DRAWINGS 

5 The file of this patent contains at least one drawing executed in color. Copies of this patent 

6 with color drawing(s) will be provided by the Patent and Trademark Office upon request and 

7 payment of the necessary fee, 

8 FIG.l illustrates the principle of relative pressure calculations for a given seismic horizon 

9 on a stacked time section. 

FIG.2 shows a flow chart of the method for calculating anomalous pressure estimates and 

ipl corrections for reflection traveltime distortions using a set of stacked seismic sections processed 

Ul2 with the aid of 3D techniques. 

ffls FIG.3 is an example of anomalous pressure map for producing sandy-argillaceous Jurassic 

1^4 formation. 

Hs FIG.4a shows an isochron map constructed using a standard technique of reflection horizon 

,£}6 tracking. 

17 FIG. 4b shows a map of traveltime corrections for variation in signal frequencies along the 

18 reflection horizon. 

19 FIG. 5 shows the final anomalous pressure map laid over isochron map and exploratory 

20 well locations. 

21 FIG. 6 depicts a three-layer model of the medium of a subterranean formation having 

22 discrete structure. 

23 FIG. 7 depicts two qualitative graphs of the parameters of fluid flow in a thin discrete 

24 layer: graph "a" depicts pressure gradients in the presence of two-axis anisotropy; and graph "b" 
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1 depicts pressure gradients in the presence of triaxial anisotropy. 

2 FIG. 8 depicts the method for determining the fluid-dynamic parameters on a seismic 

3 section. 
4 

5 DETAILED DESCRIPTION OF THE IP^NTION 

6 Taking into account the information known which is discussed in the Background section 

7 above J we concluded that each layer or stratigraphic unit of a sedimentary basin is subjected to 

8 variable load along its extent. Li one location the layer will have a high load (for example, equal to 

9 the normal lithostatic or gravitational pressure), while in the other location the load can be 
y 0 significantly lower, thus leading to the presence of a decompressed zone where reduction of general 
Ml rock pressure is associated with fracture opening. Our invention makes use of this observation 
I J2 since one may predict what will happen to the fluid occupying the pore space of a continuous rock 
W as well as interconnected system of fractures. That is, the fluid, particularly its light fractions, will 
!:14 migrate from the regions of relatively high pressure toward the decompressed zone, and if the 
y 5 closure conditions exist, accumulation and preservation of the accumulation of any type of fluid 
.^6 will take place. It should be noted that the presence of simple anticlinal structures is not necessary, 

17 and an accumulation can be formed in absolutely non-traditional conditions. The main factors in 

18 this case are the degree of discreteness of a given layer (the spatial density of fractures) and relative 

1 9 change in the general pressure. 

20 Hence, we determined that if one knows the relative distribution of the general stress within 

21 all the sedimentary sequence or within its producing interval, it is possible to establish an objective 

22 criterion for discovering subsurface fluids' accumulations independently from all known structural 

23 and lithological conditions favorable for such accumulation. Seismic data is the primary source of 

24 information for making the stress condition estimate. At the same time all existing methods and 
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1 techniques for seismic data interpretation are based on classical solutions to a wave equation for 

2 continuous and static media. In these solutions it is common to substitute stresses in initial wave 

3 equation by strains using the Hooke's law, i.e. usmg the elastic moduli. Therefore, the acoustic 

4 impedances or parameters derived from them are in practice the informative parameters for seismic 

5 data interpretation. But if the discrete dynamic model is used, the wave equation solution should be 

6 obtained for prestressed media, and the components of stress condition of real media are the target 

7 parameters. These parameters are of principal importance for exploring and discovering new oil and 

8 gas pools and water reserves as well as for planning their development. 



9 hnmediately following is a description containing the solution to a classical wave equation 

QO for a simple case of wave incidence on the interface between two half-spaces one of which is in 

;il prestress condition. We v^U now show our derivation of an expression for normal-incidence 

Yp reflection coefficient for an elastic wave: 

g3 Z-Z* + AZ'^/o) 

J4 Ar = ( 3 ) 

m z+z-^ + AZ'^/o) 

m 

where A r is reflection coefficient, Z and Z* are acoustic impedances for unstressed half-space 

'¥9 and prestressed half-space correspondingly, and co is circular frequency, 

20 AZ^ = PZ''o)/ju'' (4) 



21 where P is a value proportional to non-hydrostatic pressure in the stressed half-space, and is 

22 shear modulus in the stressed medium. 
23 

24 Reflection and refraction coefficients for plane elastic waves propagating in a prestressed 

25 medium 

26 1. Motion equations 
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1 The main law of dynamics m has the following form for a continuous elastic 

2 medium: 
3 

5 p ^= (a) 

6 dr ax J 

7 
8 

9 where stress tensor divergence in the right part of the formula is equal to the force acting on the 

10 unit volume of an elastic body subjected to internal dynamic stresses. 

1 1 The generalized Hooke's law can be written as follows: 
12 

□ l3 gui 

Tij = c ijkiS ki^c ijki (b) 

ri 17 Here c ijui is the tensor of elastic moduli. Substitution of the formula (b) into the equation (a) 
: 18 results in 

lJl9 

320 a^Ui a^ui 

P =Cijkl (c) 

^22 dXjSXk 
"23 
24 

25 For an isotropic solid body 
26 

27 c ijki= X5 ij dki + ^ ( 5 ik S ji + S ii S jk ) (d) 

28 We then assume that an elastic medium is found under the condition of static prestress, 

29 caused by external forces, defectiveness of the medium, and other reasons. We further assume 

30 that prestresses are distributed in the medium continuously. In this case elastic wave propagation 

31 in a prestressed medium can be described as superposition of small dynamic strains to, in the 

32 general case, finite static strains. The linearized theory of elastic wave propagation in a 
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prestressed medium is correct under this assumption, and consequently elastic waves can be 
described by equation (3) where the tensor of elastic moduli is given by 

c'ijki = Cijki + a^ij Sik (e) 
- where c /y k i are effective elastic moduli describing response of a prestressed elastic medium 
to small dynamic loads, and a^ij is prestress tensor. By this means we have derived linearized 
motion equations for a prestressed elastic continuous medium in the form 

P ^ =c Hki (f) 

dt dx j dx k 



2. Problem formulation. 

If one considers two elastic half-spaces having a tight planar contact, the elastic moduli 
c ijki of the lower half-space which is not subjected to prestresses differs from the elastic 
moduli of the upper half-space because of different composing materials. Besides, material of the 
upper half-space is prestressed and described by the prestress tensor 
cr ^ij. In this way the elastic moduli of the upper half-space have the form (e). We then consider 
a plane wave incident on the interface from the lower half-space. Our problem then has the 
following formulation: it is required to determine propagation direction, polarization, and 
reflection and refraction amplitudes for known propagation direction, polarization, and amplitude 
of the incident wave and specified elastic moduli of both half-spaces. 
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#2 
13 
14 
15 
16 



18 

19 




1 

2 XI 
3 

~- 4 

g 5 Boundary conditions for a tight contact require strain and stress continuity to be met at any 
iij 6 moment of time and at any point of the contact surface specified in this case by the equation /"^ 
7 x'^ = 0 , where - is the vector normal to the contact surface, and jc"^ - is coordinate in the 
^ _ 8 contact plane X2OX3. The conditions of strain and stress continuity can be written as follows 



I R T 

u i Z u i = I u i 

R T 



R T 



(g) 

(h) 



Here we use the following indices: I- for incident waves, R -for reflected waves, and 
17 T- for refracted waves. 



We then consider the following type of monochromatic waves 



(i) 



20 where m = uopi , m - is wave amplitude,/?/ - is polarization vector. 

21 The strain continuity condition (g) has the following form for the wave (i) 



E:\Files\Trans Seismic\TR0I-0Q4-PAI.doc, 12/01/99 



12 



1 

: - ^> 

6 It follows from the formula (j) that at any moment of time 

7 

^ a = to = CD (k) 

9 

10 and at any point of the contact surface 
11 

12 k-^x-" =k^ x-^= k-^"^ x-^ (1) 
13 

14 We thus derive the following formulae from the equations (1) 

15 

16 (k-^ - k'^)x-^ = 0, (k^^'-k^Jx^ =0 

t3l7 By comparison with the contact surface equation /"^ x~^ = 0, one will see that that the vectors 

•rl8 (k'^^ - k'^) and (k'^^ - k~^) are perpendicular to the contact surface, 

r 1 9 All of the wavenumber vectors lie in the incidence plane specified by the normal /"^ and 

ig20 wavenumber vector k~^ . Projections of all wavenumber vectors to the contact plane are equal to 

1=42 1 each other. We thus get 

a2 k'^sine^ =k^sm0^ =k^sin0^ 

Since the frequency content of the reflection wavelets is the same, the directions of wave 

24 propagation are determined from the following relations: 

25 sind^ /V(e^) = sine^ /V(9^) = sine ' / V (9^) 
26 

27 Now we derive simplified boundary conditions for strains from the formula (j) 

28 

29 

30 Ui + j; u i = 2^ u i (m) 

31 ^ ^ 
32 

33 Substitution of the formula (9) into the Hooke's law (b) gives the following 
34 

35 Ti / = -ic ijki'^ Ukkiexp /i(a)t -k'^x ] 
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The boundary conditions for stresses (h) now have the form 

Cijki [ k\ V* + 2■/t^ V, ; =c^,•,/2■ k\ u\ (n) 

R 

Here the elastic moduli of stressed medium c ijki are determined from the formula (e). 

To simplify further our conclusions we assume that internal stresses determined by the 
initial stress tensor cr ^/y do not violate the medium isotropy, i.e. the tensor a "ij is isotropic: 
y =± P 5ij. Here P- is the tensor diagonal or a value proportional to the pressure (sign + 
corresponds to relief, and sign - corresponds to compression) in a stressed medium. 

3. Particular case of a plane shear wave with horizontal polarization (SH-wave). 

For further conclusions we specify propagation direction and polarization of the incident 
wave. To simplify computations I select SH-wave as an incident wave. With this selection only 
two waves will appear the contact surface, i.e. reflected wave and refracted wave with the same 
polarization as that of the incident wave. So we have 
m = '^upi = ^udi3 

for each wave. Substitution of this expression into (m) results in 

0 1 OR _ 0 T 

u + U = U (O) 

Condition/; = - Sji is satisfied for the normal /"^ , and " Uk = " u Sks 
Now instead of (n) we derive 

CiJSlfk'l'u' + kWj=C^,J3lkVu'' (p) 

Since wavenumber vectors lie in the plane defined by the normal / and the 
wavenumber vector k^^ of the incident wave, they have only two components (in the axes 7 
andX2). By performing summation with respect to the mdex / we convert the expression (p) to 
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1 

2 

3 Ci}31 [kl U k 1 U ]-^Cil32 [k2 U + k 2 U J = 

4 

H" TOT + TOT 

5 ^ C ii 3 1 k I U ^ C (132 k 2 u (q) 

6 

7 Taking into account our assumption about initial stress tensor isotropy (cr%* = ±P 5 ij % we 

8 derive equation (4) into the following expression for the elastic moduli contained in equation (q) 



10 C 1131 = il Ssi + jJL ( 5 iS dll + 3 il S 13 ) 

1 1 C il32 = XS il SS2 fJ. ( S i3 Sl2 S i2 S 13 ) 
12 

Q3 It can thus be seen from the latter that only the elastic modulus is not equal to zero for an 

;H4 isotropic medium. So instead of expression (q) we derive 

^ r fl 0 I , jR 0 R ^ + jT 0 T /V 

IM C3m[ kl u + k 1 u ] = c 3131 k 1 u (r) 

. 

1 8 We now have a system of equations 

(01 OR 0 T 

^go \ u + u = u 

m i 

Vn I cnn[kW^ + yt^/V; =C3m ft""/ V (s) 

^§3 

We next introduce the reflection and refraction coefficients by definition 

25 

26 def "k* 

27 Ar 



28 
29 
30 

31 def « 

32 = 

33 

34 

35 Then we have the following solution to the equation system (s) 
36 
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1 At = Ar + 1 

2 (t) 
3 

4 C' 3131 k^l - C313I 1^1 

5 Ar = (U) 

6 C 3131 k 1 - C 3131 k 1 
1 

8 Taking into account that c 3131 = // in an isotropic medium, and using the formula (e) and 

9 a^ij = ± P S ij derive from (u) 



10 



11 



//* k^i - ^ k^i ±Pk^i 



12 Ar =Ar(P)= (v) 

13 fi k^i - fi ""k^i i±Pk^i 
14 

15 

R I I I T T T 

16 By substituting k i=-k i =k cos 0 , k i =-k cos Q we obtain 

17 

18 

19 n k^cos d - ( fi* ±P) k^ cos 0^ 

20 Ar = (w) 

21 COS 6>^+ * ^+ P) k^ cos 0^ 
22 

23 

24 After introduction of acoustic impedances of shear waves Z 3= p V3, Z*3 = p *V*3 , we obtain 
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1 J & ju (opfVs)^ 

2 kju = = = CO Zs 

3 Fi V3 
4 

5 k%* = (o Z*3 

6 Introducing notation 

7 AZ* = Pk^ 

8 we finally obtain: 
9 

10 Zs cos0^- (Zs* ±AZ*/(i) ) cose^ 

12 Zs cos + (Zs* i±A Z*/o) ) cos 0^ 

13 

=015 The formulae (u)-(x) show that the reflection coefficient and consequently the refraction 

[16 coefQcient depend on the stress condition of the other haltspace. Similar results can be obtained 

Q7 for compressional waves and for the case where both half-spaces are prestressed. Formulae for 

1 8 reflection coefficients for the said variants differ fi-om the formula (x) only by constants. 

iij9 

I Jo Inventive results of the derivation: The Method and its Application 

H\ The expression (3) is different from the classical one by an additional member {A Z* /a) 

22 which depends on the difference between stress conditions of the upper and lower half-spaces. 

23 Taking into account the fact that in a discrete medium the elastic moduli have little effect on the 

24 changes in stress condition, the expression (3) can be rewritten in the form: 

25 AZ*/oi) 

26 Ar = ( 5 ) 

27 2Z+ AZ*/(o 
28 

29 If we assume that the elastic moduli do not vary along given reflecting boundary, while pressure P 

30 varies, and the reflection coefficients in two points £ and iW of the reflecting boundary are given by 
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1 

2 R A Z^La)M(2Z-^ aZ'^m /com) 

3 = (6) 

5 

6 Taking into account the expression (4) the formula (6) can be rewritten in the following form with 

7 accuracy up to a constant: 

9 dP=—^ ; (7) 

10 Pm (A''K)'(G>Lf 

11 

12 Where: a = 2 for the media with high moduli (carbonate basins), and a = 1 for the media with 
13 

14 low moduli (terrogenous basins). 
15 

_^16 By this means the value dP is a relative estimate of pressure changeability along a seismic 

^a? reflector. For a stacked time section the amplitude of signal is proportional to the reflection 

coefficient for given reflecting boundary, and we modify expression (7) to contain analytical signal 

parameters: 

^20 

m dP^(to = (8) 

H3 (a^(ti))'' (dF^/dt)^ 

^^4 

26 where a(t) = (u(t) + v(t)f^, F(t) = arctg (v(t) /u(t)) is Hillbert transformation, and a(t) and dF/dt 

27 are instantaneous amplitudes and frequencies correspondingly. Instead of substitution of the 

28 parameters at the point M we use the following idea: amplitude parameters a(ti) and dF^/dt 

29 smoothed over the whole horizon represent a model of normal lithostatic pressure, or which is the 

30 same, "automatically" account for unknown elastic moduli of the medium. 

31 We then apply expression (8) in the method for estimating pressure gradient on a stacked time 

32 section under the following assumptions: 
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1 - aV/i) and dF^/dt are determined as average values of the instantaneous amplitudes and 

2 frequencies for the whole reflecting horizon or some its part taking into account little variability 

3 and influence of the elastic moduli of the rocks within a local part of the basin comparatively to 

4 the scale of sedimentary process, 

5 - if dP(ti) is greater than unity, this indicates a reduction of the general pressure at a given point 

6 with respect to the average (normal or gravitational) pressure, which indicates a region of 

7 unloading, 

8 - if dP(ti) is less than unity, this indicates an increase of the general pressure up to the normal 

9 (average) value or, which is less probable for discrete media, up to a value exceeding the 
Q 0 normal pressure. 

:5l FIG.l illustrates the principle of calculating parameters used in the formula (8). Time section 

^2 for the seismic line 10 is an example of standard 3D processing in the area of detailed exploration 

=^3 for oil Graph 10a shows geometry of the tracked reflecting horizon "C" coincident with the top of 

144 producing stratigraphic sequence in the terrigenous basin (sandy-argillaceous Jurassic formations). 

Q5 The fragments 16 and 18 show intervals of two traces 2 (Tm) and 5 (Tn) corresponding to said 
horizon at the points 12 and 14. Also shown here are the results of Hillbert transformation of these 

17 traces: 4 and 7 - instantaneous amplitude (A)y and 3 and 6 - instantaneous phase (F). The points / '"c 

18 mdt"c indicate position of the tracked horizon "C" after a standard picking procedure via the 

19 extrema of reflected wavelet. 

20 To obtain stable values of parameters contained in the formula (8) the following steps are 

21 performed: 

22 - find a singular point on the phase curve F closest to the point / '"^ (the curve F crosses zero 

23 at the point t"^ cr), 
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1 - find the previous singular point on the phase curve F closest to the point t*" cr (the curve 

2 F crosses value of 1 80 degrees at the point t*"^ cr), 

3 - find the following singular point on the phase curve F closest to the point t"" cr (the 

4 curve F crosses value of 1 80 degrees at the point t'"^cR/ dt), 

5 - in the time window specified as described above the normalized values of instantaneous 

6 amplitude and instantaneous fi-equency are given by 
7 

8 t'"2cR/dt 
9 

10 JJUi 

11 t"'^CR/dt 

H2 A'^c = (9) 

(t'"^CR- t"'^CR)/dt 

ills 

H7 t'"'^cR/dt 

!§8 Z(Fi-FiH)/dt 

:i9 r^cR/dt 

1=50 /'"c = (10) 

m\ (t""cR- t"'lcR)/dt 

©2 
I=t3 

■%4 * here dt is tibie trace sampling interval. 

25 The parameters A"" c and c are calculated for all points tc of the horizon "C" for a 

26 given section. As a result we obtain two sets of the normalized values of instantaneous amplitudes 

27 and frequencies and after that calculate average values of A'^c md/c using the follovwng formulae: 

28 

29 k k 

30 A\ A'c )/k ; fc-(Zfc)/k (11) 

31 i=l i=l 

32 

33 * here k is the niimber of traces in the seismic section. 
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1 Finally, use the formula (8) to determine pressure gradients (i.e. relative changeability of pressure) 

2 in each point of given horizon ^'C": 

3 dP'(t'c) = (A'c/A'cf(f''c/fcf (12) 

4 Assumedly the exponent "a" for a terrigenous basin is equal to unity. Finally we determine 

5 new values oft' c corrected for variation in the frequency spectrum of reflected wavelets along the 

6 horizon "C". As is evident from the FIG.l, the differences 

7 Atm=t'^^ cR-t"^ CR and ^ = / ''^ c/? - ^" a? for two traces 12 (/m) and 14 (w) are not equal 

8 to each other because of difference in frequency content of the same reflection in two distant points 

9 of the tracked horizon. Consequently, all the values / ' c should be corrected as follows: 
|lO t'^'^c - t'c-2(l/fc-l/fc) (13) 

J^ll It should be appreciated that the map of pressure gradients which is generated using Eqn. 

pl2 12, an example of which is shown in Fig. 3, can be of value in and of itself for identifying regions 
of differential rock pressure. Such can be used for example in the study of earthquakes, including 

u predicting areas of likely seismic activity. 
H,15 FIG.2 represents the complete chart of the proposed seismic data interpretation method 

%\6 targeted at determining relative pressure estimates and locating oil and gas fields. Preferably the 

17 input data for the interpretation method are stacked time sections obtained as a result of standard 

18 2D or 3D processing. In a first embodiment the interpretation scheme is applied to the results of 3D 

19 processing. The input data are, for example, the set of stacked time sections (gathered in in-line or 

20 on-line direction) marked in the FIG.2 by blocks 20, 40 and 42. The process, preferably, starts from 

21 the section 22. At first, a selected horizon is picked in manual or automatic mode using a generally 

22 accepted principle of tracking via the extrema of wavelets belonging to this horizon. The picked 

23 traveltimes / 'c are, preferably, saved in the matrix 34 and transmitted to the following step 24 
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1 where the Hillbert transformation is preferably performed. 

2 At the step 26, the picked traveltimes t c are advantageously used to calculate normalized 

3 values of instantaneous amplitudes and frequencies A^c and fc using the formulae (9) and (10). 

4 Then, at the step 28, we advantageously determine the average values A^c and/^c for obtained data 

5 sets. At the step 30, desired values of relative pressure changeability dP ^ are advantageously found 

6 and the obtained data set is preferably saved in the matrix 36. At the step 32 we advantageously 

7 determine corrections for signal frequency variability and correct the values f ' c preferably 

8 using the formula (13). The corrected traveltime values are preferably saved in the matrix 38. 

9 Then all the before mentioned steps are preferably applied to the section 40 and so on until 
2lO the steps have been applied to all sections. When the last section (42) has been processed, we will 
%\ have filled matrices 34, 36, and 38, which can be used for generating the following maps: initial 
l=d2 isochron map, pressure gradient map, and corrected isochron map. These maps are the basis for 
:oi3 identifying the most probable locations of oil, gas, water and other fluid accumulations. 

7;14 FIG.3 shows a relative pressure changeability map (44) for the producing horizon " C 

";15 The map has been created using the described scheme for an oil field in the terrigenous basin. The 

pJ6 regions of anomalously low relative pressure values, of which 80, 82 and 84 are representative, 

17 have been shaded red and represent a hydrocarbon indicator for fractured reservoirs in sandy- 

18 argillaceous Jurassic formations. The obtained result is then compared with the reflector geometry 

19 which is shown in FIGS. 4a and 4b. 

20 FIG. 4a shows an isochron map for given horizon (fragment 46) and FIG. 4b shows a map 

21 of traveltime corrections dt c (fragment 48). Several positive anomalies of traveltime corrections 

22 can be seen in FIG. 4b. If we apply the traveltime corrections, new anticlinal structures (namely 50, 

23 52, and 54) not previously detected on the results of standard horizon picking appear on the 

24 traveltime map. Vertical closures of the newly deUneated structures are of 20-30 m, the value of 
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1 great importance for the region of low-relief oil-bearing structures. 

2 FIG.5 shows final interpretation results: the relative pressure changes map (FIG.3) has been 

3 laid over corrected isochron map shown in the FIG.4a. The isochrons have been drawn as contours 

4 while the value of the relative pressure changeability estimates have been indicated using a color 

5 palette. The red shadmg (e.g., 58 and 68) corresponds to anomalously low values of changeability, 

6 while blue shading (e.g., 90) corresponds to normal lithostatic pressures. 

7 Analysis of the map of FIG. 5 brings us to the following conclusions: 

8 - low-pressure regions (58, 60, 62, 64, 66, 68 and 70) are localized; 

9 - these regions coincide with highs on the isochron map of FIG. 4a; 

10 - the regions 58,60, and 62 coincide with previously discovered oil fields which verifies the high 

1 1 quality and fidelity of our method for predicting the presence of fluid; 

12 - the regions 64, 66, 68 and 70 should be therefore recommended for exploratory drilling. 

13 At the final step, the depth map for producing the horizon is advantageously constructed by 

14 implementing a standard time-to-depth conversion technique to the corrected isochron map (48). 

15 This map together with pressure estimates and recommended exploratory drilling locations can be 

16 used at the prospect evaluation stage as well as at the field development planning stage, because the 

17 estimates of relative pressures in the reservoir can have the key role while selecting methods and 

1 8 technologies of oil production, like in this particular case. 

19 It is to be appreciated that the method described can, as suggested in the above paragraph, 

20 further include using classical, and future, discovered, seismic interpretation techniques in 

21 conjunction with the resulting pressure-isochron map. For example, classical seismic interpretation 

22 techniques are useful for identifying faults and traps which can cause fluids to accumulate in a 

23 subterranean reservoir. When an area of interest identified by the pressure-isochron map 

24 corresponds with the presence of a fluid trap, the likelihood that fluids are present in the identified 
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1 area is increased. 

2 The described methods are preferably performed using a computer. Thus, the invention 

3 further includes a computer configured to carry out the disclosed methods. The computer is 

4 configured to read a computer readable medium or computer memory which contains computer- 

5 executable steps for performing the steps of the method shown and described in the chart of Fig. 2. 

6 Generally, sets of computer-executable instructions for performing seismic data processing and 

7 generating resulting maps or seismic sections for interpretation are common in the industry, and in 

8 fact constitute the typical manner by which seismic methods are carried out. Examples of computer 

9 readable medium include a hard drive, a tape drive, a diskette, and ROM (read only memory) 

10 devices, such as a compact disk. Examples of computer readable memory include a hard drive and 

1 1 RAM (random access memory). As a result of performing the computer-executable steps on an 

12 appropriate input data set, the computer produces a product data set. The product data set can be 

13 outputted in the form of charts or graphs which are indicative of relative pressure changeability dP ' 

14 of a subterranean region, pressure gradients in the region, corrected traveltime values / ^ c of the 

15 region, or any of the following maps: an initial isochron map, a pressure gradient map, a corrected 

16 isochron map, or an overlay of the relative pressure changeability map and the corrected isochron 

17 map. As indicated earUer, these maps are the basis for identifying the most probable locations of 

18 oil, gas, water and other fluid accumulations. 
19 

20 Method For Predicting The Dynamic Parameters Of Fluids In A Subterranean Reservoir 

21 It is known that oil and gas pools in a sedimentary basin are formed as a corollary of three 

22 main interconnected processes: generation, migration and accumulation. Existing technologies 

23 for studying reservoirs using geological and geophysical data to identify such oil and gas pools 

24 are targeted at detection of probable locations of fluid accumulation, known as "traps". Current 
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1 methods for trap detection do not necessary lead to success in identifying such traps, however, 

2 due to many technical, geological, and other factors. However, we have developed a different 

3 approach to identifying traps having the probability of fluid presence, comprising estimation of 

4 the parameters of fluid flow or migration flow. Specifically, by determining the main parameters 

5 of the migration flow, for example its speed and direction, then the anomalous values of these 

6 parameters can correspond to the location of a trap. This method allows for the detection of traps 

7 where prior detection efforts, using attributes such as structural images of the reflecting interfaces 

8 and acoustic impedance, would be unsuccessful. The method of the present invention allows for 

9 parameters of migration flow to be determined by determining the relative estimate of the total 

10 ground pressure on the attributes of selected seismic signals. 

1 1 The invention thus further includes a method for determining the relative estimate of the 

12 total ground pressure dP on the attributes of the seismic signals related to a reflecting horizon, or 

13 any selected interval of the seismic section. In terms of geodynamics in a subterranean reservoir 

14 or basin system in question, we start by assuming that: 

15 P{'l)=Pg±dP{'l) (Eqn.l4) 

16 where: 

17 P(T) - is the total ground pressure at given depth and at given point in geological time; 

18 ^g' ^^ normal lithostatic pressure; and 

19 dP^lt) - is the supplementary pressure caused by the basin motion at the same point in 

20 geological time. 

2 1 Thus, by obtaining a relative estimate of the total ground pressure dPit) for a reflecting 

22 horizon using the seismic data, we then have a method for calculating the parameters 

23 
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1 

2 V(T),C and jj on the basis of Equation ( 1 4) using certain assumptions about the model of the 
3 

4 fracture-porous volume of the medium in the studied interval of the subsurface and the ratio of 

5 solid and liquid components of the total supplementary pressure ciP{ 7) . 

6 Recognizing the generally accepted scheme of destruction of a sedimentary basin, we can 

7 assume that the response of a discrete medium to the basement movements caused by current 

8 geodynamics of the lithosphere is formed in a slowed rate in the form of a viscous flow. 

9 Consequently, we can make a principal assumption about the property of continuity of solid and 

10 liquid components of the total pressure and further consider dP(T) as over-hydrostatic pressure 
Ol 1 which we will notify as Pd . 

:tl2 This approach allows us to solve constructively the problem of fluid migration as applied 

u l3 to the geodynamic history of a basin's development. In this case the speed and direction of a fluid 

roi4 flow V are expressed in the following form in correspondence with Equation (14): 

gi6 V(T)=-(c/ju) VdP(T) (Eqn. 15) 

y^ll where: 

19 F f - is the vector of the speed of the fluid flow; 

20 ^'is the permeability of the medium; and 

21 ^ - is the fluid viscosity. 

22 Once we have obtained a relative estimate of the total ground pressure P(T) for one of 

23 reflecting horizons using the seismic data, then we can employ a method for fmding the 

24 parameters 
25 

26 V ( T) ,c and ju using Equation (15) by making certain assumptions about the model 

E:\Files\Trans Seismic\TR0I-004-PALdoc, 12/01/99 

26 



1 of the fracture-porous volume of the medium in the studied interval of the subsurface. The 

2 derived set of the fluid-dynamic parameters of the reservoir within an oil or gas field will reflect 

3 the dynamic state of the fluid mixture at the time of determination of these parameters. 

4 Some general peculiarities of fluid migration and permeability in a discrete system are 

5 depicted in Fig. 6. Fig. 6 shows a simple and most general variant of a composition of three 

6 formation units 61, 62 and 63 of low level (B i \ and each formation unit represents a layer 

7 composed of elementary blocks. It is supposed that the geometry and orientation of these blocks in 

8 layers can be arbitrary in the plane XY. In such a scheme the total permeability c / of each layer 

9 will be determined by several main factors: 

10 The average value of the blocks' void interval in a horizontal direction ( m v ). 

1 1 The thickness of given layer ( /z ). 

12 The permeability of each block in a layer (block or so-called matrix permeability). And, 

13 The permeability of the upper and lower boundaries of a considered formation unit (boundary 

14 permeability or hydrodynamic contact along the horizontal boundaries of formation units in the 

15 context of discrete medium). This permeability is conditionally determined by the structure of a 

16 contact plane with the vertical "thickness" m h . 

17 In a generally accepted scheme of determination of true fracture permeability it is 

18 assumed that the block and fracture voids compose a common hydrodynamic system where the 

19 model of Darcy's flow is true. The permeability of a composition of layers with block structure 

20 in the normal (i.e., vertical) direction will be significantly lower than in a horizontal direction, 

21 since the value of void interval between the layers (horizontal contact of the blocks), which we 

22 notified as m a , is always considerably less than m v. Intersections of the sides of the blocks 

23 belonging to adjacent layers (in separate points such as the point S in Fig. 6) also carmot 
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significantly increase the vertical permeability of discrete rocks. Therefore, assuming some 
general direction of fluid flow exists, the speed and direction of the flow V will be different in 
each considered layer. That is, discrete systems possess an evident anisotropy of the fluid flow 
parameters, and consequently, the permeability in such systems is a tensor. 

By analogy with a generally accepted representation of flow as a conductivity, we will 
consider such a thin anisotropic layer and define the continuity equation for an incompressible 
fluid with the viscosity ju as the initial one 



vv=o 



(Eqn. 16) 



Then, let : 



c c 
V = -VP,= — 



-I 



dx, dy 



d 

J 



(Eqn. 17) 



where c is the permeability tensor in the main axes: 



c = 



-XX 
0 



0 



(Eqn. 18) 



Then fi-om the contmuity equation Eqn. (16): 



dP 



dx 



+ 



d_ 



d^P 



= c 



d 



d^P 



d 



yy 



dy^ 



= 0 



(Eqn. 19) 



Next we introduce the coordinates: ^ = 



d^P^ d^P^ 

Then: ^ H = 0 . This gives us: 

df drf 



X 



1= 



y 



'XX. 



-yy 
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In 



f 2 

X 



y 



(Eqn. 20) 



2 2 

X y 

2 The contours of equal pressure have the following form: 1 = const 

^xx ^yy 



2 '^yy 

3 and are the ellipses with the ratio of semi-axes: / = 

^xx 



4 The constant A of Eqn. 20 we obtain from the expression: 



Q=\Dndl 



(Eqn. 21) 



6 where / is the outline surrounding the region of the pressure source, and Q is flow rate of a fluid 



7 from the deformed region of small size Dn. From this we obtain the following results 



A = 



■xx^yy 



(Eqn. 22) 



9 and 



10 



In 



^xx'^yy 



2^ 



( 2 
X y 

+ — 

\^xx ^yyj 



(Eqn. 23) 



1 1 and further: 



12 



13 and 



14 



Vx = 



-xx 



dP 



d_Q 1 



H dx 2n ^c^cyy y_ 



-xx '-yy 



Vy = 



'yy^d 



Q 1 



y 



H dy 2n ^c^cyy ^ y^ 



■^xx '^yy 



(Eqn. 24) 



(Eqn. 25) 
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V X 

The vector of the flow speed V is directed radially: 77^ = — . Then we have- 

Vyy y 



And for the circle x = we have: 

Thus, it follows that the regions of higher permeability coincide with the regions of higher fluid 
speeds. But the pressure gradients are equal at ^ = 0 ; ^ = — : 



^(^ = 0) = --^^— (Eq„.28) 



and 



It IS therefore accurate to say that the pressure gradients are maximum when ^ = ^ . At 
any other (p the pressure gradients have lower values. Therefore, the circle diagram 71 for the 
pressure gradients have a form shown in Fig. 7(a). At the same time, in the presence of triaxial 
anisotropy, as shown in Fig. 6, we have a circular graph 72 for a certain plane passing through 
two main directions as depicted in Fig. 7(b). As applied to the considered problem of fluid flow 
in a discrete medium, the above description indicates that a real formation unit in the form of a 
discrete layer is characterized with the anisotropy of permeability and the ability to deviate the 
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1 fluid flow from the initial direction and to change the rate of the fluid flow. 

2 Next we consider the situation where the sources and drains of a fluid are formed in each 

3 unit volume of a discrete medium due to variation of the total voidness as a result of continuous 

4 change of the total ground pressure towards its increase (compression) or decrease 

5 (decompression). Then we have from the equation of continuity: 

d 



a 



(i^r) + V(rr) = 0. (Eqn.30) 



7 For an incompressible fluid with a constant density p it follows that 

8 ^^"""^""^ (Eqn.31) 
■0 9 where V is the speed of Darcy's flow, F is the total voidness (fracture-porous volume), and q is 
J^lO the rate of its change ( c"' ). We will also accept that the rate of the volumetric deformation 
111 1 /^/At of a certain volume of the discrete medium D is related to the rate of change of the total 
:\12 voidness solely. That is: 

% AD 

if ^^^^-'^^ 

By substituting the expression for the speed of Darcy's flow V=c/ju grad Pd in Eqn. (15), we 

15 obtain the Poisson's equation for pressure in the following form: 

17 In the case of a medium with a constant permeability cg and a homogeneous fluid with a 

1 8 viscosity juo have: 



19 APr = 



0 
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where dDo'^dxodyodzo is the volume being deformed. 

The expression for the pressure in the presence of a heterogeneous region D, in which c 
and ju are the function of coordinates, can be written as: 



V^P^^dD (Eqn.35) 



d 0 4;r^CQ// 

where r = ^ + (y^rj )j + ( z- 0k , and dD = d^djjd^. 

Equations (34) and (35) can be directly generalized for the case of several deformed and 
heterogeneous regions, and the volumes D and Do can coincide partially or completely. If we 
assume that a small deformation does not cause a considerable displacement of the permeability 
boundaries of considered volumes, there are several deformed volumes Do=ZDoj, 
j'='l,2y3....,n, and each volume Dojis characterized with its own rate of volumetric deformation 
qj . The rate of volumetric deformation has different signs in the directions of compression and 
decompression of the volumes. Accordingly it is possible to construct the sedimentary cover in 
the form of considered model of a discrete dynamic system. In a given case the problem of a 
basin's dynamics is directly associated with the problem of variation of the medium permeability 
as a time function. This eventually leads to formation of absolutely natural communicating 
pattems of the blocks: one system of blocks squeezes the fluid out, while the other system of 
blocks absorbs it. 

For regions sufficiently extensive along the axis OY in comparison with the dimensions 
along the other axes we can write the equation in two-dimensional case, similar to Equations (34) 
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1 and (35), in the form: 



13 



J 



(Eqn. 36) 



3 The integral equation for the gradient of potential can be obtained for inner points of the cross- 

4 section S using Eqn. 36 as follows: 



-1 



(Eqn. 37) 



6 where: V^=i5/a:c+ka/az, W =i d/dC+k d/dC, , dS=dldC, ,^={x-xo)\^(z-zo)k, 

7 r= (x-^\+(z-(^\i, and dSo=dxodzo. 

=C 8 By solving the integral Equation (37) we can determine the pressure from Eqn. (36) using 

9 the same matrix of inner values ^ jyP^ ■ 

Wo Now we will consider a thin horizontal layer located in a medium impermeable for a 

y 1 fluid. We will also superpose the regions of deformation and heterogeneity S and So- In this case 

\A2 Equations (36) and (37) should be written in the form 



'd 2n: 



Ins 



CJJ^ 



1 



J 



14 and 



0 



'0^ 



-1 



(Eqn. 38) 



(Eqn. 39) 



15 where Va=\ dldx+\ d/dy. VM=i d/dC+] d/dy] , dS^d^di] , and r= (x-^i+(y'7j)k . 

16 After scalarization the vector integral Equation (39) can be considered as a system of 

17 integral equations of Fredholm of the second type, if the distribution c/ju is known and the 

18 pressure gradient should be calculated, as well as the equation of Fredholm of the first type for an 
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1 unknown distribution c/^, if the distribution of the pressure gradient is specified. If the main 

2 target is evaluation of the correspondence between the pressure and the permeability, then it is 

3 proper to solve the direct problem of determining the pressure gradient on given distribution of 

4 the permeability by using Eqn. (39) as the Fredholm equation of the second type, and then to 

5 calculate the pressure using Eqn. (38). 

6 From the practical point of view the general scheme of determination of the main fluid- 

7 dynamic parameters can be presented with the aid of a diagram shown in Fig. 8. With reference 

8 to Fig. 8, the initial information used to obtain a solution to the problem is the relative 

9 distribution of over-hydrostatic pressure P d 81, which can be determined from the seismic 
^pJO attributes from the studied interval of a seismic section (2D-3D) 85 in the manner described in 
Lql 1 the first part of the above Detailed Description. Then, by solving Equation (39) with respect to an 

unknown distribution c/ ju, we can will derive c 82 in the same points and coordinates of the 

. 13 seismic section for a specified (for example, ^1). An alternative solution at this step can be 

obtained by specifying c/^ from geological data, if the location of the stratigraphic boundaries 

'M5 is known and the properties of rocks are determined on well data 83. At the final stage 84 the 

16 components Vx and Vy of the vector of fluid flow speed are derived in correspondence with the 

17 ratios firom Eqn. (26) using known distributions Pd and c/jj,. In given example the method 

18 supposes determination of the fluid-dynamic parameters on one seismic section and construction 

19 of a map in vertical plane (x - the direction of the seismic section; y - the vertical axis of time or 

20 depth). The same can be done for many seismic sections (i.e., a three-dimensional or 3D 

21 solution), and in this case horizon-oriented maps of target parameters are created for specified 

22 depths. The derived set of the fluid-dynamic parameters of the reservoir within an oil or gas field 

23 will reflect the dynamic state of fluid mixture at the time of determination of these parameters. 
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This method thus allows for the following information to be obtained: 

1. To locate the regions with highest permeability within a reservoir, or the regions with 
anomalous deviations of the vector of fluid flow speed from the dominant direction. In 
both cases such regions can testify to the presence of hydrocarbon traps. 

2. To delineate more accurately an oil or gas accumulation within the limits of given 
reservoir. 

3. To determine the discrepancy between the new and the previous outlines of the 
accumulation in the case where its spatial location was known from the results of the 
previous investigations (seismic or other methods). 

4. To optimally design the trajectory of the wellbore. This is of particular advantage for 
horizontal drilling since the highest productivity will be obtained in the situation where 
the axis of the wellbore is located at a certain angle to the direction of the maximum 
permeability ^, which is a tensor value in dynamic systems. And, 

5. To optimally plan the development scheme for given field with account for the current 
and the future dynamics of the fluid mixture. 

While the above invention has been described in language more or less specific as to 
structural and methodical features, it is to be understood, however, that the invention is not 
limited to the specific features shown and described, since the means herein disclosed comprise 
preferred forms of putting the invention into effect. The invention is, therefore, claimed in any of 
its forms or modifications within the proper scope of the appended claims appropriately 
interpreted in accordance with the doctrine of equivalents. 
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We claim: 

1. A method for deteraiining the location of the accumulation fluids in a subterranean 
formation, comprising: 

determining a first velocity vector "Vx" for migration of fluid in a region of interest in the 
subterranean formation, the first velocity vector comprising attributes of speed and direction of flow 
of fluid in a first direction in the region of interest; 

determining a second velocity vector "Vy" for migration of fluid in the region of interest, the 
second velocity vector comprising attributes of speed and direction of flow of fluid in a second 
direction in the region of interest; 

extrapolating the velocity vectors to identify the fluid accumulation location; and 

wherein the first and second velocity vectors are primarily functions of supplementary 
pressure "dP" in the region of interest, the permeability "c" of the region of interest, and the 
viscosity "u" of the fluid in the region of interest. 
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2. The method of claim 1 wherein the supplementary pressure is determined by identifying 
pressure gradients within the region, said region being characterized by a seismic image, said 
seismic image comprising a stacked time section representing horizons within said region, 
comprising: 

a) picking a first selected horizon from said seismic image; 

b) calculating a set of instantaneous amplitudes and frequencies for said first selected 
horizon; 

c) determining the average amplitude and frequency of said set of instantaneous 
amplitudes and frequencies; 

d) identifying pressure gradients associated with said instantaneous amplitudes and 
frequencies to generate a pressure gradient map, said pressure gradients corresponding to points 
at which said instantaneous amplitudes and frequencies vary from said average amplitude and 
frequency, wherein points at which said instantaneous amplitudes and frequencies are less than 
said average amplitude and frequency correspond to locations of relatively low pressure. 

3. The method of claim 2 wherein said first selected horizon has associated traveltimes, and 
wherein said instantaneous amplitudes and frequencies are calculated by the Hillbert 
transformation using said traveltimes. 

4. The method of claim 3 wherein said pressure gradient associated with said traveltime 

dP^ (t^ c) is calculated using the formula 
dP' (t' c)-=(A' c/a\/ (f^c// cf . 
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5. The method of claim 1 wherein the velocity vectors are related to the permeability "c", the 
viscosity "u", and the of supplementary pressure in the region of interest dP by the equation 

V(T) = -(c/^) VdP(T) 



The method of claim 1 wherein the permeability "c" is calculated for selected values of 



the permeability "u" using the equation 
1 r-^O r 



w^here Va=\ dldx+] d/dy, Vm=i d/dC+} d/dr] , dS=d^dr] , r is a scalar = (x- 
^i+(y-Tj)k, and n are incremental lengths in the respective directions x any y, 



A = 



and Q is flow rate of the fluid in a portion of the region of interest. 



7. The method of claim 1 wherein the permeability "c" and the viscosity "u" are obtained from 
geological data in the region of interest. 
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1 8. The method of claim 1 wherein the first velocity vector "Vx" is calculated using the equation 



1 Q \c^r 

y ( __r^ 5 and the second velocity vector "Vy" is calculated using the 



-12. 



equation V^^^_^ = ~ 2^ y — ' wherein Q is flow rate of the fluid in a portion 



4 of the region of interest. 

5 
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1 ABSTRACT 

2 A method for determining the location of the accumulation fluids in a subterranean 

3 formation. The method includes the steps of determining a first velocity vector "Vx" for 

4 migration of fluid in a region of interest in the subterranean formation. The first velocity vector 

5 includes attributes of speed and direction of flow of fluid in a first direction in the region of 

6 interest. The method further includes determining a second velocity vector "Vy" for migration of 

7 fluid in the region of interest. The second velocity vector includes attributes of speed and 

8 direction of flow of fluid in a second direction in the region of interest. The velocity vectors are 

9 then extrapolated to identify the fluid accumulation location. The first and second velocity 

10 vectors are primarily functions of supplementary pressure "dP" in the region of interest, the 

1 1 permeability "c" of the region of interest, and the viscosity "u" of the fluid in the region of 

12 interest. 

13 The supplementary pressure can be determined by identifying pressure gradients within the 

14 region, the region being characterized by a seismic image of a stacked time section representing 

15 horizons within the region. The permeability of the media within the region, and the viscosity of 

16 the fluid within the region, can either be determined mathematically or from geological data. 
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Fig. 8. 



